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ABSTRACT 


ANALYSIS OF PULTRUSION PROCESSING FOR LONG FIBER 
REINFORCED THERMOPLASTIC COMPOSITE SYSTEM 


Pultrusion is one of the composite processing technology, commonly recognized as 
a simple and cost-effective means for the manufacturing of fiber-reinforced, resin matrix 
composite parts with different regular geometries. Previously, because the majority of 
the pultruded composite parts were made of thermosetting resin matrix, emphasis of the 
analysis on the process has been on the conservation of energy from various sources, 
such as heat conduction and the curing kinetics of the resin system. Analysis on the 
flow aspect of the process was almost absent in the literature for thermosetting process. 
With the increasing uses of thermoplastic materials, it is desirable to obtain the detailed 
velocity and pressure profiles inside the pultrusion die. Using a modified Darcy’s law for 
flow through porous media, closed form analytical solutions for the velocity and pressure 
distributions inside the pultrusion die are obtained for the first time. This enables us 
to estimate the magnitude of viscous dissipation and its effects on the pultruded parts. 
Pulling forces refined in the pultrusion processing are also analyzed. The analytical model 
derived in this study can be used to advance our knowledge and control of the pultrusion 
process for fiber reinforced thermoplastic composite parts. 
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Chapter 1 
INTRODUCTION 

1.1 Preview and Literature Survey 

The demands for stronger, tougher, higher temperature resistant thermoplastic com- 
posites has increased dramatically in recent years. Thus, the pultrusion technology faces 
many new challenges. Since 1957, over 130 papers dealing with subjects related to the 
pultrusion process, such as product and tooling designs, raw materials, machinery, prod- 
uct markets, etc. have appeared. However, few investigations dealt with in-depth analysis 
of the flow field inside the pultrusion die. Because of the lack of basic understanding of 
the stress and strain fields generated by the pultruded laminate part inside the heated die, 
most pultrusion processing was still operated on a trial and error basis. Such a practice 
leads to frequent fiber breakages and process down-times. 

The previous works focused on the physical and chemical properties of thermosetting 
resin systems. For this type of resin system, prediction of the temperature inside the 
pultrusion die requires accurate reaction kinetics and material properties. The equations 
of reaction kinetics for a variety of chemical systems are often nonlinear. When coupled 
to the heat transfer to and from the die wall, the resulted equation for conservation of 
energy is difficult to compute. 

Price [l] 1 was the first to use a heat-transfer model for pultrusion analysis. Two 
limiting cases were examined: an isothermal case with a uniform die wall temperature 
and adiabatic case where heat conduction was negligible. The model used first order 


The number; in brackets indicate reference 
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kinetics for epoxy resins. No experimental data was provided for the evaluation of 
modeling results, however. 

Tulig [2] used finite elements to model pultrusion cure of epoxy resins in round and 
irregular die shapes. Boundary conditions was specified to simulate both the heat input 
from the die heater and heat losses due to convection with air. Tulig’s work is the only 
published model to date which has been successfully verified with experimental data for 
epoxy resin. 

Han, et al. [3] used an autocatalytic model for unsaturated polyesters and epoxies, 
and allow density, thermal conductivity and heat capacity to change with degree of cure. 
No experimental data was presented, however. 

Ma, et at. [4] published a model similar to Han, but, for the first time, axial 
conduction along the pultrusion die was included in the calculations. No evidence was 
given, however, to show that axial conduction was significant. 

Batch and Macosko [5] used a mechanistic kinetic model for polyesters which 
included provisions for diffusion-limited chain propagations. Also included in their 
analysis were models for pressure and pulling force predictions. 

The analysis of pultrusion processing currently undertaken and presented below 
differs from the works mentioned above in two major ways: 

(1) This research deals with thermoplastic based composite systems. There is no 
reaction heat evolved from the processing. In addition to the conduction heat from the 
die wall, the viscous heating generated in the high shearing zone between the pultrudate 
and the die is also included in the analysis. 

(2) Unlike the conventional assumption, that the pultrusion is equivalent to a plug 
resin flow reactor, used in literature, detail analyses of the flow and pressure fields inside 
the pultrusion die are conducted. This analysis also enable one to estimate the extent 
of viscous heating mentioned in (1) which was unable to be performed by all previous 
workers. 
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1.2 Outline 

The flow analysis for both the entrance section and the main pultrusion die section 
is presented in Chap. 2. The detailed pressure distribution is recognized as one of the 
most significant factors that determine the effect of process. It is also necessary to know 
that the relationship among the back flow pressure p b incoming volume fraction V 0 and 

permeability k. 

Chapter 3 covers the motion of the resin and fiber composite. The complete force 
distributions are given in detail by different mathematical formulas. As a matter of fact, 
the viscous drag F v , frictions Ff and collimation force F c , are affected by the processing 
parameters, geometry parameters and material parameters. 

The viscous heating is what we most expect to know in the thermoplastic pultrusion 
thermal analysis, distinguishing from the reaction-dominated thermosetting processing. 
Consequently, the conservation of the equation is derived in Chap. 4. 

The finite element analysis is used to analyze the thermal aspect of the pultrusion 
processing, which is presented in the Chap. 5. The conclusion and the appendix are 
complemented to discussion in the previous chapters. 



Chapter 2 

FLOW ANALYSIS OF PULTRUSION PROCESS 
2.1 Statement of the Problem 

The pultrusion process line consists of several processing steps. They are fiber 
impregnation, preheating zone, pultruding die and part cutoffs. Among them, the 
pultruding die is the focal point of this study. 

The pultrusion process of continuous fiber reinforced polymeric resin matrix laminate 
through a cylindrical die is illustrated in Fig. 2.1. The die consists of two sections: a 
short tapered section with length oL near the entrance and a main pultrusion die section 
with length L. A Cylindrical Coordinate Systemis is selected with origin fixed in the 
inlet of the main pultrusion die section. The composite laminate is pulled by a force F 
entering the die at z = -oL with a constant speed Vf. The contraction ratio of the tapered 
entrance section is 1/e. The main pultrusion die section is mildly tapered with inlet radius 
R at z = 0 and exit radius AR at z = L. The value for A = 0.98 is used throughout this 
investigation. It is desired to obtain the velocity and pressure profiles inside the entire 
pultrusion die, -aL s z ss L, such that the pulling forces required for the processing of 
thermoplastic composite laminate can be better estimated analytically. 


4 
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2.2 Flow Analysis in the Entrance Section 


2.2.1 Assumption 

The flow inside the entrance tapered section of the die can be analyzed as a flow 
through porous fiber bundle media with the following assumptions: 

1. Portions of the incoming fluid volume which are pultruded through the die, i.e., 
excluding the backflow volume, travel with the same line speed of the fiber bundles 
throughout the entrance section. This is a plug flow assumption. 

2. Backflow volume is estimated solely from the volume contraction of the flow by 
the tapered die geometry. 

3. One dimensional ( in -z direction ) Darcy’s law is applied, i.e., v z = v z (z) only 
for the backflows. Because the entrance section is short, secondary flows, v r (z,r), due to 
the tapered geometry, is negligible. This allows to obtain a simple analytical estimate 
of the backflow effect. 

4. The process is isothermal and in a steady state. There is a uniform viscosity 
without chemical reactions. 

5. Following dimensionless terms are defined as 

D, = Df/R, V b , = V bl /V,,k j = k,/R 2 ,P b = PbR 2 /(/iLv,), 

2 = z/I,, P T = Ptke, P T = P T R 2 /(k c /iLv f ), P en = P b + P T 

2.2.2 Back Flow Pressure Distribution Pb(z) 

The cross section area of the tapered die at z can be obtained from die geometry as 


and 


r(z) = R — z 


(c — 1)R 
aL 


( 2 . 1 ) 


A(z) = 7T 


R — z 


(e-l)R 

OrL 


( 2 . 2 ) 
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where —ah s: z s: 0, and e > 1. One dimensional Darcy’s law in the — z direction for 
the backflow volume is given by 


- < Vbz(z) > = 


k z (z) dP b (z) 


(2.3) 


fi dz 

where the averaged velocity is — < v/, z >= ~Qh(z)/A(z). The backflow volume Q b (z) 
is estimated at any location z ( -cvL szsO) based upon assumptions 1 and 2 as 


Qb(z) = Q(z) - Q(0) 

=v f [l - Vf(z)]A(z) -v f [l - Vf(0)]A(0) 

(2.4) 


where V f (z) is the fiber volume fraction and v f is the pultnision line speed. It is noted 
that the backflow volume vanishes at z = 0, which is the inlet of the main pultrusion 
die section. Thus, we have 


< v bz (z) >= - 


Qb(z) 

A(z) 


= v f { [i - V f (0)] 


A(0) 

A(z) 


(l-Vf(z)]} 


(2.5) 


Equation (2.5) was also independently derived by Batch [6], Because of the tapered 
geometry, the axial permeability, k z (z) [16], due to the change of the fiber volume 
fraction, Vf(z) is given as 


_ D; (1 - V,( z )) 3 

l( J 16Sz Vf(z) 2 


with 


Vf(z) = Vo 


A(-a)l 


( 2 . 6 ) 


(2.7) 


L A ( z ) J 

where Vo =Vf(-a) is the incoming volume fraction of the fiber bundles. The quantity Df 
is the fiber filament diameter and S z = 0.7 is a characteristic parameter of the graphite 
fiber bundles investigated by Gutowaski [13]. Values of S z have also been reported to 
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be 0.68 by Lam and Kardos [14], 0.48 by Williams et al. [15]. Generally, higher values 
of S z are observed for close packed arrays of fiber. 

Using the dimensionless terms defined before, and neglecting the underscore bars, 
the relation for the backflow pressure is expressed as 

Z 

P h (z) = -J ~ > dz (-«<*<°) < 2 8 > 

— or 

where 


< v bz (z) > = 



(2.9) 


k(7) = _gL a-vtw )! 

M ’ 16Sz Vf (z) 2 


( 2 . 10 ) 


V f (z) - V 0 


' A(-a) ' 
. A(z) . 


V 0 



( 2 . 11 ) 


A(z) 



(e — 1)R 

z — — 


with the boundary condition 


( 2 . 12 ) 


z — —o , Pb(— or) = 0 (213) 

It is noted that the dimensionless pressure distribution, P b (z), is determined by the 
permeability, k 7 (z), only for a given tapered die geometry at the entrance. In order 
to solve Eq. (2.8), values of V 0 , e, a, D f and S z have to be first specified. 

With R = .125 inch, D f = .0005 inch and S 7 = 1.0, the behaviors of P b (z) and 
k 2 (z) in the entrance section of the pultrusion die as a function of process variables are 
tabulated in Table 2.1. It is seen that the permeability, k 7 (0), at z = 0 is a function of 
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the contraction ratio, Eq. (2.11). For a given die geometry defined by a set of a and e, 
values of k z (0) decrease with increasing V 0 . The backflow pressure, Pb(0), at z = 0 is 
directly proportional to the tapered die length, a, and increases with increasing e and V 0 . 

An example of change in fiber volume fraction, Vf(z), along entrance section of the 
pultrusion die is shown in Fig. 2.2. Values of Vf is noted to increase rapidly in this 
typical die geometry. 
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Table 2.1 Changes of Pb(0) and k 7 (0) in the Entrance Section of 
the Pultrusion Die as a Function of Process Variables Indicated* 


€=1.5 a=0.1 

c=1.5 a=0.2 

e=1.5 a=0.3 

V 0 

MO) 

Pb(0) 

V 0 

k 2 (0) 

Pb(0) 

V 0 

MO) ; 

Pb(0) 

.05 

.00005 

193 

.05 

.00005 

387 

.05 

.00005 

580 

.1 

9x10 6 

993 

.1 

9x10 6 

1986 

.1 

9x10 6 

2979 

.2 

Rigl 

7202 

.2 


14405 

.2 

an 

21607 

e=2.5 a=0. 1 

e=2.5 a=0.2 

c=2.5 a=0.3 

V 0 

M0) 

Pb(0) 

V 0 

kz(0) 

Pb(«) 

V 0 

kz(0) 

Pb(0) 

.02 

.000042 

143 

.02 

.000042 

288 

.02 

.000042 

432 

.05 

.000003 

1259 

.05 

.000003 

2520 

.05 

.000003 

3779 

.1 

mm 

11173 

.1 


22347 

.1 


33521 


*R=0.125 in., D f =0.0005 in. and S z =1.0. 








































































Fiber Volume Fraction V f (z) 
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2.2.3 Thermally Induced Pressurization Distribution Pj(z) 

For a thermoplastic laminate pultruded through the die, the resin matrix temperature 
rises continuously, from the ambient temperature T am |, to the pultrusion temperature T pu |, 
due to the heat conduction from the die wall. Thermal expansion of the resin system 
beyond the confinement of the pultrusion die will give rise to a pressurization effect. The 
amount of the pressure rise due to thermal expansion is related to resin compressibility. 
By assuming that resin temperature is T = T(r, z), with constant coefficient of thermal 
expaasion a v and compressibility k c , Batch [6] derived the following expression for the 
thermally induced pressurization effect: 

P T (z)= — [T(z)-T amb ] (2.14) 

K c 

where the averaged resin temperature is given by 

h 

T(z) = p- J T(z,r)rdr (2.15) 

o 

If we further assume that the resin temperature is uniform over any cross section of the 
die along z, where — oL <. z <. L, and reaches die temperature linearly from To(-oL) = 
T am h to T(0) = T pu | at the end of tapered die entrance section (z=0), the dimensionless 
thermally induced pressure Pj ’(z), as defined above, can be estimated by 

P T (z) = a v (T p „,-T» mb )(l + j) (-a<z<0) (2.16) 

In order to compare the magnitude of dimensionless pressure contributions on an equal 
basis, a new dimensionless term P T is required as 

P T — P T R 2 /(/c c ^Lvf) 

Consequently, Eq. (2.16) becomes: 
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+ (in, 

A new dimensionless pressure term is also defined as P en (z) = P b (z) + Et( z ) 
sums up the pressure generated in the entrance tapered die section in the pultrusion 
process. Using the set of typical values for the pultrusion parameters in Table 2.2, 
magnitudes of P^(z) and P-p(z) are plotted in Fig. 2.3. It is seen that the thermally 
induced pressurization is a dominant factor for the pressure build-up in the tapered 
section near the entrance of the pultrusion die. The basestone of making a successful 
mathematical model is the selection of the parameters which are refined from the 
experimental data. The type of materials can be categorized as materials parameters, 
geometry parameters, and processing parameters with given dimensionless or dimensional 
values. 
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Table 2.2 Typical Values of Material, Die Geometry 
and Processing Variables in the Pultrusion Processing 


[ Material parameters ~T 

Geometry parameters 


n 

1.5 

a 

0.1 

h( oo) 

0.87 

A 

0.98 

V f ( oo) 

0.74 

e 

1.5 

/T 

0.8 

R 

6.35x 10‘ 3 m 

■ 

Kfenec 

2.72x 10 7 

L 

0.3048m 

V 0 

0.25 

Processing parameter 

s 

.x 

D { 

0.002 

TpuJ-Tamb : 

200° C 

S z ~ 1 

1.0 

Vf 

5.08x 10 3 m/sec 

n 

0.14 PaSec 



r 

a v 

5x 10' 5 °C l 




5x 10' 10 m 2 N 1 






Pressure Distribution P(z)xlO 
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2.3 Flow Analysis in the Main Pulrusion Die Section 


2.3.1 Assumption 


The flow inside the main pultrusion die section can be analyzed as a flow through 
porous fiber bundle media with the following assumptions: 

1. The inlet of the main pultrusion section is the end of the tapered section near 
the die entrance. Following the assumption 1 specified for the analysis of the tapered 
section, we have a flat inlet velocity profile with a magnitude of v z (0,r) = v f , the same 
as the pultrusion line speed. 

2. The conservation equations are satisfied analytically. There is no transverse or 
longitudinal variations of volume fraction for the resin ( Vf = const ). 

3. The contraction ratio of the die, A=0.98, close to unity, is used in this study. This 
is a necessary geometry to obtain a complete analytical closed form solution for the flow 
field inside the die due to the mildly tapered geometry which is, however, negligibly 
small. 

4. Unlike epoxy and polyester matrix resin, the impregnated resin matrix used in 
the current investigation is thermoplastic in nature with Newtonian properties. There is 
no chemical reactions occurring during the processing. Pressurizations due to thermal 
expansion, vaporization and shrinkage of the curing resin inside the die (commonly 
occurred for the thermosetting resin system) are absent. 

5. The process is isothermal and in a steady state. Heat transfer is neglected, i.e., there 
is no temperature gradient inside the die and the resin matrix has a uniform viscosity [7]. 


6. The one directional permeability assumption in the main die section is character- 
ized by the following relation (see also Fig. 2.1) 


M z ) = 


Df (l-Vf(z)) 3 

16Sz V f (z) 2 


(2.18) 


where 
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h(z) = Az 4- 1 — z 

m i 

M z ) [Az + 1 — z ] 2 


(2.19) 


( 2 . 20 ) 


V f (z) = V,(0) 

= v,(0) 


A(0) 


A(z)J 


[Az -f 1 — zj 


( 2 . 21 ) 


Mz) 

M°) 

By knowing values of Vf(0)(or k 7 (0)), one can calculate k z (z) at any z along the main 
die by Eqs. (2.21) and (2.22). The permeability for the unidirectional fiber structure 
exhibits anisotropic properties. For simplicity, the permeability of isotropic material is 
considered as constant, coincidently with the assumption 2. 

7. Following dimensionless terms are defined: 
r = r/R, z = z/L, k = k/R 2 , v = v/vf, P = PR 2 /(jxLvf), 
f(r) = f(r)/(v f /L),h = h/R = 1 - z(l - A) 

2.3.2 Governing Equations and Boundary Conditions 

An empirically modified Darcy’s law suggested by Brinkman [8, 9] together with the 
condition of incompressibility for the flow through porous media is: 


v,(0) 


1 - V,(z) 

V,(z). 


l-v,(0). 


( 2 . 22 ) 


v - k V 2 v = (VP) 


(2.23) 


V • v = 0 


(2.24) 
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When the additional term with Laplace operator is dropped, Eq. (2.23) resembles the 
Darcy’s law. An obvious difficulty of using the ordinary Darcy’s law in calculations of 
flow problem is that the viscous shearing stresses acting on a volume element of fluid have 
been neglected; only the damping force of the porous mass has been retained. Hence the 
ordinary Darcy’s law cannot be used as such in detailed flow analysis. Equation (2.23) 
with the inclusion of viscous stress tensor term has the advantage of approximating the 
ordinary Darcy’s law for low values of k. When the values of k are high, it approximates 
the Navier Stokes equation for the viscous flow in empty space. The additional term 
with Laplace operator in Eq. (2.23) was intended to account for distortion of the velocity 
profiles near die walls. Unlike the ordinary Darcy’s law, when Eqs. (2.23) and (2.24) 
are applied to flow through a porous medium in a tube, the result can be simplified to 
the Hagen-Poiseuille law when k approaches infinite [10]. 

Because of the assumption of small taper for the main pultrusion die, i.e., A = 0.98, 
we may assume that v 7 = v 7 (z,r), v r = v r (r), v 0 = 0 and P = P(z). Consequently, Eqs. 
(2.23) and (2.24) become 


19 / dvA 

r 8t V dr ) 


k r 


1 d 
rd7 ( 


\ 3 2 v,l 

k z dP 

(2.25) 

+ 
CD ' 

N 

^ t 

II 

1 

/t dz 

H| I* 

II 

o 


(2.26) 

8v. 

+ dz 0 


(2.27) 


Note that the secondary flow, v r (r), is assumed to be a function of r only. From Eq. 
(2.27), it is obtained that v 7 (z,r) = zf(r) + Vf, and d 7 v z /dz 2 = 0 and, therefore, Eq. 
(2.25) becomes 


df 2 ldf 1 dP Vf 

dr 2 ^ r dr k Zfi • dz ^ zk 


(dimensional) 


(2.28) 
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Equation (2.28) is a dimensional form where the permeability in z-direction is denoted 
by k = k z = k z (0), which was obtainable from flow analysis conducted earlier for the 
tapered section near the die entrance. Using the definitions of dimensionless terms, the 
dimensionless form of Eq. (2.28) with the underscore bars neglected is expressed as 

iLI _l 1^1 _ If = - f— 4--^ =Co (dimensionless) (2.29) 

dr 2 rdr k z\dz^k7 

where C 0 is a constant. Equation (2.29) has the form of the modified Bessel’s equation 
for f(r) with n = 0, a = KT ,/2 (see Appendix A). 

2.3.3 The General Solution in the Dimensionless Form 

A general solution of Eq. (2.29) is readily available as [11] 


fM = C, !„(-)=) + c ’Ko(^) -kCo (2 ' 30) 

The dimensionless boundary conditions are specified as 
r = 0, df/dr = 0, 
r = h, f = f(h) = -1/z 

Thus, the complete ( dimensionless) solution of Eq. (2.29) is found be 


f(r) = I(v z (z,r) - 1) = 

z 


kC 0 

z 


Io (jg) 

L io (^k) J 


-kC 0 


(2.31) 


Consequently, the velocity distribution inside the main pultrusion die section is expressed 
as 


v z(z,0 = [1 - zkCo] 



which satisfies the following boundary conditions 


(2.32) 


r — 0, dv z jdr — 0, 
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r — h, v z (z, h) = 0 

Note that v z (0,0) = l-l/lo(l/k ,/2 ) approximates 1 ( or in dimensional form v z (0) approx- 
imates Vf ) from Eq. (2.32). Because of the small taper of the die, the secondary flow 
v r (r) is expected to be negligibly insignificant. 

From Eq. (2.29), we have 


dP 

dz 



The value of constant Co is found to be 


(2.33) 


Co =2 


l ~ p ' n(0) 


(2.34) 


Solution of Eq. (2.33) is a parabolic function, and it is expressed as 


POO 


1 

k 


Pen(O) 


rz + P en (0) 
k 


(2.35) 


This satisfies the boundary conditions: P(z) = P en (0) at z = 0 and P(z) = 0 at z = 1, 
where P en (0) is the pressure at the end of the tapered section near the die entrance, as 
calculated before. The pressure distribution is inversely proportional to the permeability 
k. Substituting Eq. (2.34) into Eq. (2.32), we have 


v 2 (z,r) = {1 


2[1 — kP en (0)z]} 



(2.36) 


where the dimensionless h = l-z(l-A). 
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2.3.4 Discussion 


Figure 2.4 shows the pultrusion die pressure distribution, Eq. (2.35), for various 
values of P en (0) indicated with k = 9.19x10^. From Eqs. (2.33) and(2.34), we have 


dP 

dz 


= 2 



z 


1 

k 


(2.37) 


It is noted that at z = 0, dP/dz is negative, and other cases are discussed as follows: 


Condition (i): P en (0) = 0, then 0 < dP/dz = 1/k at z = 1. In this case, a minimum 
pressure occurs within the pultrusion die at z m j n = 1/2. 

Condition (ii): dP/dz > 0 or 0 < P cn (0) < 1 /2k at z = 1. In this case, the minimum 
pressure occurs within the pultrusion die with 


z t (2.38) 

I "' n 2[1 - kP en (0)] 

where 1/2 > z m j n > 1.0. Such a behavior of conditions (i) and (ii) was not observed 
experimentally for the pultrusion process. 

In order to describe a pressure distribution with physical meaning, it is required to 
present other cases decribed below. 

Condition (iii): dP/dz ;> 0 or l/2k 5 P en (0) at z = 1 . In this case, the minimum pressure 
occurs at (when P en (0) = l/2k) or beyond ( when P en (0) > l/2k ) the die exit (z = 1.0). 
Condition (iv): For r = 0 and any given z = z k along the die, Eq. (2.36) gives 


v*(z k> 0)={l-2[l-kPen(0)]z k }[l- 



(2.39) 


It is seen that in order to have V 7 (z k ) :> V z (0, 0) for any z k * 0, one must have P en (0) * 
1/k. Consequently, we note that P en (0) > 1/k is a necessary condition to obtain pressure 
and velocity distributions in a pultrusion die with physical significance. 
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From the set of typical values for the pultrusion process given in Table 2.2, the 
pressure and velocity profiles can be calculated. The change of permeability along the 
pultrusion die is shown in Fig. 2.5. Values of k decrease dramatically in the entrance 
section of the die as a result of fast rise in fiber volume fraction, Vf(z), shown in Fig. 
2.2. It remains relatively unchanged within the main die section. The pressure profile 
in main die section is shown in Fig. 2.4. The pressure builts up quickly in the tapered 
section of the die within short length (a = 0.1), then decreases monotonically, with 
increasing negative slopes, and eventually vanishes at the die exit shown in Fig. 2.6. In 
the present example, it is noted that P en (0)=1.8xl0 7 and k=3.39xl0~ 7 , and this results 
into P en (0)xk=6.1 which is greater than 1. Velocity distributions in three cross sections 
(z = 0.2, 0.5, 0.8) along the main die are shown in Figs. (2. 7-2.9). The velocity profiles 
closely resemble the plug flows, i.e., the profiles are rather flat except in the narrow 
regimes near the die wall. It is noted that the velocities at the fiat portions of the profiles 
are all higher than the fiber pulling speed, Vf=5.08xl0~ 3 m/sec. 

Substituting Eq. (2.36) into Eq. (2.27), the velocity v r (r) at any given location z 
along the main section of the pultrusion die can be obtained by integration as 


v r (r) = [1 - kP etI (0)]r - 2Vk[l - kP en (0)]{ 

. It 


(1 - A){1 - 2(1 - kPen(0)]z k }{ 


'■(ve) 

>•(&) 


}+ 


Io 


(2.40) 


which satisfies the boundary condition of v r (0) along the pultrusion die. The velocity 
v r (r) at any given z = z^ and h^ = l-z^t-A) can be shown to be negligibly small when 
compared with v 7 (zk, r). 




Fig. 2.4 Pressure Distribution in the Main Die; a=0.1, e-1.5, V o -0.1 and k 9.19x10 
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Chapter 3 

ANALYSIS OF PULLING FORCES 
3.1 Calculation of the Pulling Force 

Forces generated in a die are functions of stock velocity, reinforcing configuration, 
resin system and pultruding conditions, materials, etc. They are summarized below 
according to the acting effects: 

(a) Frictional forces ( fiber friction against the die wall ). 

(b) Viscous forces ( shear viscous flow in a thin layer ). 

(c) Collimation forces ( backflow drag resistance on fibers, fiber compact ). 

(d) Temperture-induced forces ( increasing viscosity and resin thermal expansion ). 
Following dimensionless force terms are defined : 

F v = F v /(//Lv f ),Fj = F}/(//Lv f ),F I f I = F}V(/^Lv f ), F I { 1 = F{ n /(R- 2 LK FENE c) 

P = PR 2 /(^Lv f ),P (e = Pf e /(KFF,NECR)>KFENEC = K F ENEcR 2 /(^ v f) 

E ” 1 = ET^FENEC ’Efe = £-f e K. E ENEC >E C — Pc/(R LKfenec) 

l c c = F'K fenec ,Ei = Fc/(/rLv f ) , F k = F c k /(^Lv f ) 
h(z) = h(z)/R , h(— a) = h(— or)/R = e , h(oo) = h(co)/R 
where , ^F 1 } 1 ^ and (F£) t are written without tilde in Table 3.1. 
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3.1.1 Viscous Drag Force Contribution, F v 

We first consider the viscous drag contribution in the entire pultrusion die — aL s 
z s L. This is expressed as 


Li 

F v = J r„] r=h 2jrh dz 


— e»L 


“ / ('"SL 2rhdz+ 

~al 

L 

I l —ft — -] 27rh dz (dimensional) 

J \ dr / r=h 


(3.1) 


The first integration of Eq. (3.1) between the length (-aL « z s L) is equal to zero 
because of the “plug flow” assumption made in the entrance tapered section. This 
may not be serious assumption. Using the dimensionless terms defined above without 
the underscore bar, and noting the dimensionless relationships of the main die section 
geometry h = 1 — z(l — A),dh = —(1 — A)dz, we have 



(3.2) 


Equation (3.2) is the contribution of viscous drag force from the main die section to the 
total pultrusion force. It is noted that F v is a function of permeability k and kP en (0) only 
for a given die geometry and pultruded prepreg system; the results are shown in Fig. 3.1. 
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3.1.2 Friction Force Contribution, Ff 

The friction force is originated solely from the relative movement between fiber 
bundles and the die wall, and is the product of the frictional coefficient (/{f) and the 
normal forces acting on the wall surfaces. The normal forces acting on the wall surfaces 
have four sources of contribution: (I) the hydrostatic pressure, (II) the normal stress 
generated by the flow inside the pultrusion die, (III) the fiber compaction forces due to 
the contraction geometry of the die, and (IV) the vaporization and shrinkage forces due 
to chemical curing reactions of resin system. In the present investigation we deal with 
the thermoplastic materials, effect due to the reactions do not exist in our consideration. 
The other three sources of contribution to the friction force are presented below. 


3. 1.2.1 Hydrostatic Pressure Contribution, Ff* From the information on area, normal 
force, and frictional coefficient, /q, the relation for Ff 1 is expressed as 


L 

F{ = 2 7r//f j P(z) h(z) dz (dimensional) (3-3) 

-«L 

Since a “plug flow” assumption with a uniform velocity v z (z,r) = vj was made in the 
tapered die section, there is no pressure gradient along length (-«L ^ z ^ 0). Therefore 
Eq. (3.3) is simplified as 


L 


F f = 



(dimensional) 


(3-4) 


o 

With the dimensionless forms for all parameters difined earlier, the dimensionless ex- 
pression for Eq. (3.4) without the underscore bar is as follows: 


l 

Ff = 2jr/(f J f ff)p(z) h(z) dz (3.5) 

0 

where h = 1 — z(l — A) , 0 < z < 1 and 
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P(z) = 


i- p ™(0) 


-ft) 


z + P e „(0) 


(3.6) 


Substituting Eq. (3.6) into Eq. (3.5) and performing the integration, we obtain a closed 
form solution as 


F i = <s)(f){- L ^ + kP '”<°> 



(3.7) 


It is noted that Ff 1 is a function of k and kP cn (0) only. The solution given by Eq. (3.7) is 
illustrated in Fig. 3.2. It is seen that Ff’ decrease continuously with increasing values of k. 
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3. 1.2.2 Normal Stress Contribution, Ff 11 For a given frictional coefficient, //f, the 
relation for Ff 11 is given by 

L 

F} 1 = 2nm J {r rr } r=h h ( z ) dz (dimensional) (3.8) 

-a L 

where r rr = -2/*(dv r /dr) for the Newtonian fluid is the normal stress in r-direction. 
From the continuity equation, Eq. (2.27), we have r rr = 2/z(v r /r + dv z /dz). Again the 
“the plug flow” assumption with constant v z (z,r) and negligible v r is used in the flow 
analysis of the tapered die section. Thus, Eq. (3.8) becomes 

L 

F} 1 = = 27T//,f J (r rr } r=h h(z)dz (dimensional) (3.9) 

o 

The dimensionless form of Eq. (3.9) without the underscore bar is given as 


F S' = 2 '"-/{ 2 [T + (r)^]L Mz)dh 




+ 1 LJ dz 




hdh 


(3.10) 


In the main pultrusion die section (fi s z 5 1 or A s h ^ 1), the secondary flow, v r (r), 
at the wall (r = h) is neglected. Consequently, we have 


Ff 1 = 


27T 


'»/ Ki)SL>* 


where h = 1 - z(l - A), 0 < z < 1 or A < h < 1 and 


(3.11) 


5V Z A x?V,i . 

-Jr=h = “(I - A ) Jr=h 


dz 


= — 1={(1 - A) - 2[1 - kPen(0)](l - *»)} 

Vk 


■ (^r) 

W*)J 


(3.12) 
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A combination of Eqs. (3.11) and (3.12) results in 


F ' I 211 - kP '" (0)l(1 - ' h » 


L Io (^r)J 


hdh (3.13) 


It is noted that Ff H (k) is a function of k and kP en (0) only. The results of this equation 
are shown in Fig. 3.3. 


3.1.2.3 Fiber Compaction Contribution, Ff llf A composite laminate pultruded through 
the die with a tapered entrance section is expected to experience continuous compaction 
in the lateral r-direction. The laminate could be viewed as a fiber bundle with multi- 
filament layers in which exists numerous points of contact among the individual filament 
layers. Any degree of compaction applied to the fiber bundle will raise the fiber volume 
fraction, and increase the fiber elastic forces. It was observed experimentally [12-15] 
that forces required to compact the fiber bundles increased dramatically when the fiber 
volume fraction, Vf, approached a limiting value of Considering the fiber bundles 
as a whole an elastic spring, we proposed a phenomenological model for the fiber elastic 
force, F e , based upon the Finitely Extendable Nonlinear Elastic (FENE) spring concept: 


F e (t) = Kfene{ 


h(0) - h(t) 


fi_ 

Vh(0)-h(oo)/ 


} 


(3.14) 


,h(0)-h(oo), 

where F e (t) is the elastic force at h(t), and h(t) is the thickness of the fiber mat at t 
with A(0) denoting the initial thickness with no compaction occurred. The difference 

A A 

MO) — M°°) represents the maximum compaction achievable for a given fiber mat. This 
model has three adjustable parameters: n, h(co) and spring constant Keene- A fiber mat 
with this force law will behave as a linear (Hookean) spring for small compaction, but will 
get stiffer and stiffer (nonlinear behavior) as the compaction increased. Furthermore, the 
laminate cannot be compacted beyond A(oo) (or in an other word, exceed A(0) — A(oo)), 
because infinitely large compaction force will then be required according to Eq. (3.14). 
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In the pultrusion process presently studied, total fiber volume remains constant 

A 

throughout the pultrusion die of Fig. 2.1. The quantity h(z) in Eq. (3.14) will be 
the radius of the pultrusion die at z with h( 0) denoting the radius of die entrance (ft(0) 
= h(-o) = e) where no compaction occurs. Then for a laminate rod pultruded through 
the tapered die, we have 

h(-o) — h(z) 

_ ( h(-g)-&(») 

Vh(-«)-h(oo) 

where Pf e (z) is the pressure acting on a characteristic surface of the die wall at z, and 
Kfenec = KFENE/(tmit area) in the units of (force per(length) 3 ) is a material property 
determined experimentally for a given random fiber mat or aligned unidirctional prepreg 
system. The results of Eq. (3.15) are shown in Fig. 3.4. Using the dimensionless terms 
defined earlier, we have the dimensionless form of Eq. (3.15), with the underscore bars 
neglected, as 


Pfe(z) = KfENEc{- 


yj^} (dimensional) (3.15) 


with 


PfeW - { 


h(-ct) — h(z) 


( M-a) 

- 




n} 




(3.16) 


h(z) = 1 — z 




-a < z < 0 


(3.17) 


h(z) = l -z(l-A), o<z<l (3.18) 

There are no physical meanings attached to the parameters K pen EC an< ! n - However, 
it is reported that this model can describe the deformations of both random fiber mats and 
aligned prepreg fibers rather well with reasonable values of Vo = Vy(— or), Voo = Vy(oo). 
It is seen that this model is simple but adequate for meeting the objective of present 
study, and is used in the following discussion. 
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For a frictional coefficient, /if, we have fiber compaction to the total pulling force 
in a dimensional form as 


h 

F} 11 = 2 7r /if J Pfe(z) h(z) dz (dimensional) 

— oL 


(3.19) 


The dimensionless form of Eq. (3.19) with underscore bar neglected is 


F|" = (F[") T + (F|") 
1 

J P fe (z)h(z) 


= 2 7T // f 


M 

dz 


u i 

= 2irfi( J Pfe(z) h(z) cos <f> dz + 27r/if J Pf e (z) h(z) dz (3.20) 


In order to compare the magnitudes of various force contributions, additional dimen- 
sionless force, pressure and elastic spring constant are defined as: , jP , K-eenec > 

respectively. Then, neglecting the underscore bar, the new dimensionless form of Eq. 
(3.20) becomes: 


u 1 

Ff 11 = 27r/if{ J P fe(z) cos <f>dz T J Pf e (z) h(z)dz} (3.2 1 ) 


or 


Fj n = 2n, 


l 

oKfeNEc{^-^y^ [{ | j ^ v , n }hcos<?i>dh} 

e f 1 “ V«-K(oo)J 


+ 


(-rb)/ 


<7 ‘ h x.h lhdh 


1 


(rrpby)] 


(3.22) 
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Equation (3.22) is obtained by means of Eqs. (3.17) and (3.18) and <f> = tan 1 [(e - l)/o]- 
Because of the tapered geometry, the fiber volume fraction Vf(z) is expected to increase 
along the die. As shown before for the cylindrical die (-a ^ z ^ 0), we have 


Vf(z) = Vo 


' A(-a) ' 
. A(z) . 


= Vo 


b(~ ° f ) 

. M z ) . 


(3.23) 


or 


h(-a) /VfM ( 3.24) 

h(oo) V v ° 

For an aligned prepreg system, n=1.5, V/( oo) = 0.74 and K F ENE = K FEN EC - 
50 psi, were found for the FENE model (Eq.(3.14)) to describe elastic forces reasonably 
well. For R = 0.125 in, V 0 = 0.1 and h(-cr) = e = 1.5 from Table 2.2, values of K FE NEC 
= 400 lb in ~ 3 and h( oo) = 0.55 can be calculated. 

Total contribution of the frictional forces to the pulling force is: Ff = Fj + F^ + Fj 11 
given by Eqs. (3.7), (3.13) and (3.22). The comparison among frictional forces is shown 

in Fig. 3.5. 



Pressure P fe (z) lb /in 





Fiber Compaction, F f III 
Normal Stress, Fj 
Hydro. Pressure, F^ 


Permeability, k 

a Sources of Contribution for Frictional 
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3.1.3 Collimation Force Contribution, F c 

Collimation forces are related to the geometry of the tapered die section near entrance. 
There are three different sources of contribution: (i) the bulk compaction force, (ii) resin 
backflow, and (iii) thermal expansion. They are discussed as follows. 

3. 1.3.1 Back Compaction Force, F c c As discussed before, the laminate pultruded 
through the die experiences continuous contraction in the lateral r-direction due to the 
tapered geometry near the die entrance (-qL s z s 0). Elastic force of the fiber bundles 
increases fiber volume fraction. This elastic force acts perpendicularly to the surfaces of 
the tapered die wall. We have for the bulk compaction force F c c : 

o 

F£ = (F^) T = 2x J Pf e (z) h(z) sin <^dz (dimensional) (3.25) 

— »L 

The dimensionless form of Eq. (3.25) is 

l 

K = 2,7 (- “ t) / Pfc < h > h sin<^dh (3.26) 

f 

with <t> = tan -1 ((e — I )/c] and 


Pf.W = 



t—h 

(r^y) 


} 


(3.27) 


where h(z) = 1 — z(— ■), —a < z < 0. 

In order to compare magnitudes of various force contributions, an additional dimen- 
sionless force term F_ c c is used. Then, with the underscore bar neglected, a dimensionless 
form of Eq. (3.26) becomes: 


F" = 2zr 


h)l 


Pf e (h)h sin <f> dh 


(3.28) 


or 
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F c c = 2xKfenec(-^) Jiy 7 Z rr }h9m ^ dh 
v / i i 1 - 7=wj>\ 

where values of n, a, e and h(oo) are given in Table 2.2. 


(3.29) 


3. 1.3.2 Backflow Force Contribution, F c b The dimensionless pressure profile P b (z) in 
the tapered section of the die generated by the resin backflow is given by Eq. (2.8). 
Thus, we have 


o 



(3.30) 


Using the relations for F b c and P c as given before, the dimensionl ess form of Eq 
(3.30) without the underscore bar is expressed as 


f c= 2 */ 

— Of 

where <f> = tan -1 [(e — l)/**]- 


sin <f> dz 


(3.31) 


3. 1.3. 3 Thermal Expansion Force Contribution, F c k 1 he relation for F c is given by 


F ' = ( F ') t + ( F ')m 

o ) 

= 2k f P T (z)h(z)sin^dz + 27r/t f I P T (z)h(z) cos <^dz (dimensional) (3.32) 

J n 


The dimensionless form of Eq. (3.32) is 


° /L \ } 

pk _ 2 t r j Pt( z )M z ) s ’ n \^R J P T( z ) h ( z ) COS 
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o 

“KiO / p!r(z) [ 1_ ( 1 T L ) J ] sin#dz 

— a 

1 

+27 J Pt( z )[^ ~ — ^) z l cos < f >< ^ 2 (3.33) 

o 

where <f> is as defined in Eq. (3.31). The dimensionless form of Pt( z ) is given by 


P tW = ££< T H - T„b) (l + £) (-« < * S 0) (3.34) 

Substitute Eq. (3.34) into Eq. (3.33), an analytical solution in the closed form is obtained 
as 


F: = 2tt| 


a a(e — 1) 

2 + 6 


kc)(/.v,) (TH T * mb) 

[sin <j> + fif cos <f>] + ^/i/(l + A)| (3.35) 


The thermal force contribution F c k can be estimated by the values of pultrusion parameters 
given in Table 2.2. The comparison among the collimation forces is shown in Fig. 3.6. 
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Symbolized Values of Force 

Fig. 3.6 Comparison among the Collimation Forces; or=0.1, e-1.5, 

V o =0.25, AT=200°C, D f =0.002, S z =l, a v =5xl(T 5 °CT 1 , *c=5xl0 10 m 2 N \ 
/i=0.14P a S cc , R=9.525x 10~ 3 m, L=0.203m, Vf =6.35x lCTWsec, k=9.19xl(T 6 , 
P en (0)=1.798x 10 5 , ^tf=0.8, Kfenfx= 2 - 72x1 ° 7 ’ h(oo)-0.87 and V f (oo) 
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3.2 Discussion 

Total pulling force for the pultrusion process is the sum of contributions from viscous 
drag, F v , frictions Ff, and the collimation force F c . Thus, F iota \ is given by 


Ftotal = F v + Ff + F c 
= n + F/ + F/' + (F/") r +(F/") M 
+F' + f* + (F*) T + {F*) u 


(3.36) 


It is noted that different sources of contribution to the total dimensionless pulling force 
are functions of various combinations of dimensionless processing, geometry and material 
parameters. Sets of geometric parameters which affect various force contribution sources 
are all different. Fora given die geometry and prepreg system [17], F v , Fj , Fj 1 ,and F h c 
are dependent on the permeability, k, while Fj 11 , F c c and F* are functions of the pulling 
speed, v f . The terms F 1 ^ and are written without tilde in Table 3.2. The comparison 
among the different sources of pultrusion process force is shown in Fig. 3.7. 



9 


Table 3.2 Forces with Explicit Pultrusion Process Parameters 






Sources of Pultrusion Process Force 
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Viscous Drag (3.3x10 ) 


Friction Force (1.08x10 ) 



8 


Collimation Force (2.43x10 ) 


II « 1 


i 

1 

m 

m 

m 

W 

mKKm 

rri i 

I i 

"7 

Y===^==-=^==^ ^ T 


0 5 10 15 20 Zb 

Symbolized Values of Force 

Fig 3 7 Comparison among .he Main Pulsion Process Forces; 0=0.1. <=1-5, 
V o =0.25, AT=200°C, D ( =0.002. 8,-1. Ov=5xl<T> »C', rc^xlO-'W, 
II-0-14P.S.C, R=9.525xl<T’m, L=0.203m. Vf=6.35xKT 3 m*ec, k=9.1»xl<r* 

, „ 7 72x 10 7 n=l 5 h(ex>)=0.87 and V f (oo)-0.74. 

Pen (0) = 1-798xl0 5 , /jf=0.8 , K F enec =: 2.72 x 10 l.*, ) 




Chapter 4 

VISCOUS HEATING 


The generation of heat due to the action of viscous dissipation can lead to significant 
temperature variations across the shear fields in any of viscometer configurations. For 
Newtonion fluid with temperature independent properties, with an isothermal wall, the 
temperature profile will be shown in the next chapter. This case is thought of as a 
simplest case for the study of the effect of viscous heating with constant viscosity. 


4.1 Calculation of Viscous Heating 

The viscous heating is used in high viscous flows. It is predicted by an empirical 
equation, which is a simplified form of the conservation equation of energy. For 
nondimensionalizing the energy equation, additional dimensionless terms are introduced 
as follows: 

T 0 =T/T 0 , B r (Brinkman number)=^ f 2 /(K c T 0 ), R e (Reynolds number)=pv f R//i, 
Pr(Prandtl number)=/zC p /K c . 

4.1.1 Governing Equation 

The energy equation in terms of T 0 and C p is given as 


pc p ^ = -(V • q) - (r • Vv) 

For constant p, K c , and C p , Eq. (4.1) is expressed as 


(dimensional) 


(4.1) 


DT 1 

/>c r -gp = K e v ! T + 


(4.2) 
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where Kc is the heat conductivity defined as heat flow across a unit area when temperature 
gradient is unity (Btu/(tLT)). The symbol C p is the heat capacity at constant pressure per 
unit mass (Btu/(MT)), and 7 = Vu + (Vu) 7 " . 

By noting that vr=vr(r) and T=T(z,r). Eq. (4.2) is expressed in a dimensionless 
form as 



Further, it is assumed that heat conduction in the longitudinal (fiber) direction is negli- 
gible, the velocity gradient (ov r /ar)»dv 7 /dz), and the secondary flow, v r (r), inside the 


die is negligibly small. Thus, Eq. (4.3) reduces to 

d 2 T IdT n _ /R\ dT ( dv z \ 2 _ 

dr 2 + r Or e?r ( L J V * dz + Br \ dr ) 

Equation (4.4) is solved for (dimensionless) temperature decrease due 


(4.4) 


to hpiihno 


4.1.2 Boundary Conditions 

For steady state flow condition (assumption) and T=T(r,z), the following boundary 
conditions are specified: 

At r = 0, dT/dr — 0 (dimensional) 

At r = h, T = To (dimensional) 

The dimensionless form of the boundary conditions are given as 

At r = 0, dT/dr = 0 
At r = h, T = 1 

4.2 Discussion 

This problem is solved by the finite element method. The procedures and results 
are presented in Chap. 5. 



Chapter 5 

FINITE-ELEMENT METHOD THERMAL 
ANALYSIS FOR THE PULTRUSION PROCESSING 

The finite element method is selected as a numerical approach to carry out the 
solution of the physical problem. In this method, mathematical formulations are displayed 
explicitly in a sequence of representations of the integral equations with continuous 
boundary conditions. These are different from the finite difference method which is 
applied in solving the problems where mathematical formulations are given by differential 
equations. In the thermal analysis for the pultrusion processing, the finite difference 
method (which is suited to solve fluid dynamics problems in rectangular boxes) is not 
applicable in general. In this study, the two-dimensional heat transfer equation with 
cylindrical coordinates is solved by a finite element method [18]. 

5.1 Introduction 

The analysis of thermoplastic pultrusion includes the derivation of equation for the 
velocity distribution, resin pressure, pulling force and heat transfer in the processing. 
It now remains investigation of the temperature profile in the system. The temperature 
profile can be obtained by conducting experiments. It can also be obtained by solving 
the energy equation using the Galerkin Finite Element Method. The assumption made 
in isothermal heating condition for the process leads to the boundary conditions that are 
easier to define. Also, the symmetric die design shape provides convenience for the 
temperature evaluation by the Finite Element Method. For simplicity, the assumptions 
of low viscous heating and a Newtonian fluid is made for the following analyses. The 
detailed formulations and numerical examples for the problem is given next in order to 
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demonstrate the applicability of the finite element methods for heat transfer analysis of 
the pultrusion processing. 


5.2 Governing Differential Equation 

The governing differential equations for thermal energy transfer analysis associated 
with the pultrusion processing are summarized as follows: 


^Y^ + ^ + i^_^m^ +flr ^y =0 ( 5.,) 

\Lj dz 2 dr 2 r dr \L J dz \ dr ) 

where h = 1 - z{l - A). The velocity profile v 7 (z,r) is defined along the flow direction, 

the z-direction, as 

, _ 

A ✓ . \ 


„ z (z,r) = {l-2[l-fcP en (0)]z} 


(?*) J 


(5.2) 


where 1 G is the modified Bessel function of zero order of the first kind. Velocity gradient 
along r-direction is given as 


d‘4 z ' r> = -~L(1 - 2[1 - M>„„(0)] 2 } 

dr y/k 


'(^r) 


(5.3) 


where li is the modified Bessel function of order one of the first kind. 

The boundary conditions for the low viscous heating of Newtonian fluids pultrusion 


are given as 

7’(0, r) = 1, at z=0 
T(l,r) = 1, at z=l 
T(z, h) = 1, at r=h 
|I( 2 ,0) = 0,atr=0 

where z=0 and z=l represent the locations at the inlet and the outlet of the flow processing, 
respectively; and r=0 and r=h represent the locations at the center line and the surface 
of the flow domain, respectively. 
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5.3 Finite Element Modeling 

The following discussion focuses on formulating the finite element modeling. The 
fluid domain is first divided into eight-nodes quadrilateral elements. In the finite element 
method, the field variable which is temperature in this case , is interpolated as the product 
of shape functions Nj and nodal degrees of freedom, Tj, as 


8 

T (e) = N i T i e) = NTT ( 5 ‘ 4 ) 

i=l 

where Nj, Tj are functions of r and z. 


The Galerkin Finite Element method invokes the condition that the weighted residuals 


of the governing differential equations are zero, i.e., 


/*■■'(? 

n* 


2 d 2 T d 2 T IdT _ _ 

+ -7T2+-—- R * P ' 


dz 2 


dr 


dr 


>©■*♦»(*)’ 


* = 1 , 2 , 


)dn = o, 

• • -8 (5.5) 


This leads to a matrix equation in terms of nodal temperature, Tj. Upon substituting Eq. 
(5.4) into Eq. (5.5) and integrating the resultant equation by parks, the final equation 


is obtained as 



dNi 

~dz 



T RcPr 



<m(Tj) 


n i = 1,2, • • -8 (5.6) 


The area integals defined in Eq. (5.6) can be numerically evaluated by using the technique 
of 3x3 Gaussian quadrature. Specifically, one can define the components of the local 
stiffness matrices and the force vector as 



n 


dN it ,dNj 

! 7 {r) - sT 


+ 



(5.7) 
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k}- 

'j 


= J {*’ ■ * ■ (T) ■ r ■ Ni v -ih} 

n e 


dfl 


(5.8) 


2 

• r • B r • dfl (5.9) 

O' 

where i=l 8, j=l 8. 

The technique of 3x3 Gaussian quadratures maps physical domain fi e of Fig. 5.1 to a 
square bounded by -Iss^l and -l^t^l. In this case, the shape functions are mapped to: 

Ni(s, t) = .25(1 - a)( 1 - *)(— 1 -s-t) 

N 2 (s, t) = .25(1 + s)(l - <)(— 1 +s-t) 

Nj(s, t) = .25(1 + s)(l + <)(-l + s + t) 

Ni(s,t) = .25(1 -s)(l +f)(-l -s+t) 



N$(s, t) = .25(1 - s x a)(l -t) 

N 6 {s,t) = .25(1 -f s)(l -t xt) 

N 7 (s,t) = .25(1 - 5 x a)(l + 0 
N s {$,t) = .25(1 -s)(l -t x t) 

where s, t are dimensionless coordinates after mapping transformation from r, z. Next, 
the integration of equations Eqs. (5. 7-5.9) is performed to obtain the local stiffness 
matrices as 


11 ^ 

K a — J J D(s,t)®(s,t)T)^(s,t)dsdt 

-l-l 

1 1 „ 

K b = / / N(s, <)V T N(.s,<)L(s,<)D T (5, t)dsdt 

-l -l 

and for the force vector as 

l l 

rT/ 


F = 


J J N(s,<)N T (5,<)|J| dsdt 

-l -l 


H 
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The D matrix appearing in and Kb is given as 


D = 


3N 

ds 


dN 


dt 


(1-0(23+0 

4 

(l-5)(5+20 I 

4 

(1-0(23-0 

4 

(l+a)(— j+2t) 
4 

(1 + 0(23+0 
4 

(l+a)(a+2t) 

4 

(1+0(25-0 

4 

(1 -*)( 3+2<) 

4 

-3(1-1) 

-(!-’) 

2 

tl=?l 

2 

-<(1 + 5) 

— 5(1 -f t) 

(!-’) 

2 

L 2 

-f(l - s) . 


(5.10) 


The matrix given by Eq. (5.10) represents the gradients of shape functions. The Jacobian 
matrix J is defined as 


J = 


f Td N 

_ T3N i 

« 7TT 

r t J7 

r TdN 

r T8N 

t In 

r t ~5T J 


(5.11) 


where the partial derivatives of the shape functions have been selected in terms of s, t 
coordinates as 


oz os 


(5.12) 


and 


3N T 

“dT 



(5.13) 


where Ji and J 2 are, respectively, the first and second rows of inverse of the Jacobian 
matrix J. The velocity vector V is a function of constant v z which can be written 
approximately as 


v z « N T v z 

The velocity vector V is then defined as 
V = R e ■ P r • v z • r • R/L 

The same relation can be applied to represent dv z /dr. As a result, one can define matrix 
H appearing in F as 
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H = (r • B T {d^l dr) 2 ) 

The matrix L is defined as L = Ji|J|, where |J| is determinant of the Jacobian matrix. 
Moreover, the 0 matrix appearing in K a is given as 

0 = J|J| = [jfj, ^ • r + Jjj 2 • r] | J | (5.14) 


Finally, using 3x3 Gaussian points, the integrals for K a , Kb and F can be simplified as 


K a 


3 3 

E E WiWjDufSJijDj 

i=l j=l 


(5.15) 


K b 


3 3 

£E W i W i N U VTN U L u D U 


i=l j=l 


(5.16) 


and 


E E WjWjNyNy |J|ij 


H 


i=l j=i 


(5.17) 


where subscripts i and j refer to the Gaussian points in the three-point Gaussian quadra- 
tures along s axis and t axis. The weighing coefficients Wj= W 2 = W 3 = 5/9 and Sj=- 
.774597, 0, .774597 for i=l to 3, and tj=-.774597, 0, .774597 for j=l to 3. 

After assembling the above local matrices, one obtains a global matrix equation: 


(K a + K b )T = F (5.18) 

where K a and are associated with heat conduction and heat convection , respectively. 
Standard process such as LU factorization can be used to solve the above equations for 
nodal values of temperature. 


5.4 Numerical Example 

To verify the finite element procedure proposed in the last section, a simple heat 
transfer problem is presented here. The domain of the problem is a rectangular con- 
figuration shown in Fig. 5.1, which is discretized into 12 elements and 51 nodes with 
specified boundary conditions. 




Fig. 5.1 Boundary Conditions for Numerical Example. 
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The heat transfer problem is governed by the following differential equation: 

(?) 


2 d 2 T d 2 T IdT . ,dT , r , , ft . • , a /c 1Q x 

-^T + -^T + ~!r- a r » z ) + /( r > 2 ) = 0 ( dimensionless ) (5.19) 

oz 2 or 2 r or oz 


where a(r, z) and f(r, z) are given by 
a(r, z) = R e P , 


/ d\ /o(4r) 

V ( T ) (> - 2 I> - - 7T\. 

W Mtf) 


(5.20) 


and 

f(r,z) = -^^ 4(2r 2 + l)-8(2z 2 + 3) + a(r,z)4z(2r 2 + l) (5.21) 

Equation (5.19) is very similar to Eq. (5.1) with a(r, z) and f(r, z) related to v z and 
dv z /dr, respectively. The exact solution of Eqs. (5.19—5.21) is 


T(r, z) = (2z 2 + 3)(2r 2 + l) 


(5.22) 


if the boundary conditions of the problem are specified as: 

z = 0, T(r) = 3(2r 2 + l) 
z= 1, T(r) = 5(2r 2 + 1) 
r = 1, T(z) = 3(2z 2 + 3) 
r = 0, dT /dr = 0 

Based upon the finite element formulation in Sec. 5.3, the resultant numerical temperature 
distribution is listed in Table 5.1. 
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Table 5.1 Comparison between Numerical and Actual Values 


Number 


1 


2 


3 


Exact 

Temperture 


3.000 


3.369 


4.470 


6.308 


8.880 


8.940 


9.000 


9.036 


8.916 


4.494 


3.020 


3.080 


3.454 


4.577 


6.449 


9.069 


9.130 


9.191 


9.644 


9.548 


4.813 


3.245 


3.500 


3.920 


5.182 


7.282 


Numerical 

Temperture 


3.000 


3.369 


4.470 


6.300 


8.880 


8.940 


9.000 


9.036 


8.895 


4.519 


3.586 


3.584 


3.422 


4.588 


6.499 


8.995 


8.846 


9.191 


9.644 


9.507 


4.833 


3.891 


3.835 


3.924 


5.188 


7.282 


Element 

Exact 

Numerical 

Number 

Temperture 

Temperture 

27 

10.226 

10.208 

28 

10.292 

10.278 

29 

10.361 

10.360 

30 

11.336 

11.340 

31 

11.185 

11.164 

32 

5.680 

5.694 

33 

3.845 

4.259 

34 

4.280 

4.631 

35 

4.788 

4.763 

36 

6.310 

6.307 

37 

8.854 

8.872 

38 

12.406 

12.364 

39 

12.484 

12.383 

40 

12.568 

12.570 

41 

13.530 

13.530 

42 

13.340 

13.332 

43 

6.802 

6.810 

44 

4.620 

4.853 

45 

5.000 

5.000 

46 

5.590 

5.590 

47 

7.350 

7.350 

48 

10.293 

10.290 

49 

14.417 

14.410 

50 

14.506 

14.510 

51 

14.604 

14.604 
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The column of the exact temperature consists of data obtained from Eq. (5.22). In 
the table, it is observed that the maximum difference between the exact value and the 
finite element solution is about 20% appearing at the fourth boundary. The temperatures 
obtained by the finite element method in other areas are very close to the exact solution. 

5.5 Numerical Results for Pultrusion Processing 

The finite element method developed above is applied for thermal analysis of a 
pultrusion process problem. Due to the symmetry of the problem geometry, only the 
upper half plane of the flow field is considered for finite element analysis. The finite 
element mesh is shown in Fig. 5.2, which consists of 12 elements and 51 nodes. 

Heat transfer analysis of this pultrusion process is performed and the results are 
reported in this section. Table 5.2 documents the numerical results obtained by the finite 
element method given in Sec. 5.3. 

In order to evaluate the effect of viscosity in temperature distribution, heat transfer 
analyses are performed for different Pj(z ) and B r which are the processing parameters 
related to the viscosity. These parameters are expressed as 
Pj {z) = (a v R 2 / KcfiLvf) (T pu i - T am j)( 1 + z/ot) 

and 

B T = H/{KcT 0 ) 

The numerical study shows that the final temperature distribution is uniform every- 
where for the parameters values described in Table 5.2 in which B r varied from 10 4 
to 10 4 . The factors that affect the temperature distribution in the process include higher 
curing point and superior toughness. However, these factors were not considered in the 
present case. As a result, the temperature is uniformly distributed across the die section. 
In their experimental study, Larock and Hahn [3] also observed the similar phenomenon. 
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Table 5.2 Key Processing Parameters and Calculation of Temperature Values 


viscousity 

Brinkman 

KPen(O) 

Pen(O) 

Permeability 

Temp. 

(aO 

No. (Br) 



(K) 

(T) 

0.14 

3 

1.5 

1.5xl0 6 

10 6 

1 . 






Chapter 6 
CONCLUSIONS 


Pultrusion processing for long fiber reinforced thermoplastic composite is analyzed 
theoretically. The pultrusion die consists of two sectioas: short tapered section near the 
entrance followed by a main pultrusion die with near constant diameter. The pressure 
distribution in the entrance section is analyzed, which takes into account the contribution 
from back flow and thermally induced pressurization effects. It is found that the latter 
exhibits predominant effect on the pressure built up in the entrance die section. 

Flow analysis in the main pultrusion die section is accomplished by using a modified 
Darcy’s law for flow through a porous media. The modified form incorporates viscous 
stress term which is used to account for distribution of velocity profile near the die 
wall. Closed form solutions for the velocity and pressure distribution are obtained for 
the first time. It is found that velocity profile, v z (z, r), is a function of kp en (0) where k 
is the dimensionless permeability and p en (0) is the dimensionless pressure at the end of 
entrance section of the pultrusion die. It is noted p en (0)2:l/k is a necessary condition to 
obtain pressure and velocity distribution in a pultrusion die with physical significance. 
In practice, this implies that in order to achieve a laminate flow with the pultrusion die, 
the pressure built up p cn (0) has to be greater than the inverse of permeability, 1/k. The 
velocity profiles at various die sections, z, are found to remain flat until r/R is greater 
than 0.98. 

Pulling forces encountered in the pultrusion of long fiber reinforced thermoplastic 
composite are also analyzed. Contributions from various sources, namely frictional, 
viscous, collimation are considered. It is found that viscous drag contribution to the 
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pulling forces is negligible, and collimation force contribution dominates the total pulling 
forces. 

For talcing account of viscous heating, an additional energy equation is formulated. 
Numerical analysis using finite element method is conducted to obtain temperature 
and velocity profiles inside the pultrusion die. Within the range of Brinkman number 
KT^BfidO 4 investigated, the viscous heating is found to be negligible, resulting in a 
uniform temperature distribution along the cross section of the pultrusion die. 
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APPENDIX 


MODIFIED BESSEL EQUATION 


A modified Bessel equation is 

+ xa£ - (a 2 + = o 

A complete solution of this equation is 
y = cjl n (ax) + c 2 K n (ox) for any n. 


where Iq( 0) = l,li(0) = 0,K o (0) = Ki(0) = oo. Forn-0,1,2 an integer, we have 


I n (x) = I_ n (x) = + 


+ 


2 n n! 1 2 2 l!(n + l) 2 4 2!(n + l)(n + 2) 

x 6 

2 6 3!(n + l)(n + 2)(n +3) + 


(A.l) 


Following theorems can be established 

^[x n I n (x)| = x n I n — i(x), ^[x“ n I n (x)] = x- n I n+ i(x) 
^[x n K n (x)] = — x n K n _i(x), ^[x“ n K n (x)] = — x _n K n +i(x) 
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